function th = fdr_thresholds(plist, levels, lnmeth, fdrfac)
% fdr_thresholds  - compute probability thresholds for q(FDR) levels
%
% FORMAT:       th = fdr_thresholds(plist, levels [, lnmeth [, fdrfac]])
%
% Input fields:
%
%       plist       list of p-values to be thresholded
%       levels      q(FDR) levels
%       lnmeth      method to use, default: 0
%                   0 :  c(V) = 1
%                   1 :  c(V) = ln(V) + E
%                   2 :  c(V) = 1, use last crossing
%                   3 :  c(V) = ln(V) + E, use last crossing
%       fdrfac      FDR factor (default: number of values)
%
% Output fields:
%
%       th          p-thresholds for q(FDR)

% Version:  v0.7f
% Build:    8110521
% Date:     Nov-05 2008, 9:00 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 2 || ...
   (~isa(plist, 'double') && ...
    ~isa(plist, 'single')) || ...
    numel(plist) ~= length(plist) || ...
    any(plist < 0) || ...
   (~isa(levels, 'double') && ...
    ~isa(levels, 'single')) || ...
    isempty(levels) || ...
    numel(levels) ~= length(levels) || ...
    any(levels < 0 | levels >= 1)
    error( ...
        'BVQXtools:BadArgument', ...
        'Bad or missing first argument.' ...
    );
end

% sort plist
splist = sort(plist);
levels = levels(:)';
th = zsz(levels);
nlist = numel(splist);
if nargin < 4 || ...
   ~isa(fdrfac, 'double') || ...
    numel(fdrfac) ~= 1 || ...
    isinf(fdrfac) || ...
    isnan(fdrfac) || ...
    fdrfac < 1
    fdrfac = nlist;
end

% divide with FDR factor list for fast find access
fsplist = splist(:)' ./ ((1/fdrfac):((1-1/fdrfac)/(nlist-1)):1);
lastcross = false;
if nargin > 2 && ...
    isa(lnmeth, 'double') && ...
    numel(lnmeth) == 1 && ...
    any((0:3) == lnmeth)
    if lnmeth > 1
        lastcross = true;
    end
    if lnmeth == 1 || ...
        lnmeth == 3
        lnmeth = true;
        th = [th(:), th(:)];
    else
        lnmeth = false;
        th = th(:);
    end
else
    lnmeth = false;
end

% find thresholds
for thc = 1:numel(levels)
    if lastcross
        thi = find(fsplist <= levels(thc));
        if isempty(thi)
            th(thc) = splist(1) / 2.001;
        else
            th(thc) = splist(thi(end));
        end
    else
        thi = find(fsplist > levels(thc));
        if isempty(thi)
            th(thc) = splist(end);
        elseif thi(1) > 1
            th(thc) = splist(thi(1) - 1);
        else
            th(thc) = splist(1) / 2.001;
        end
    end
    if lnmeth
        if lastcross
            thi = find(fsplist <= (levels(thc) / (log(nlist) + 0.5772)));
            if isempty(thi)
                th(thc, 2) = splist(1) / 2.001;
            else
                th(thc, 2) = splist(thi(end));
            end
        else
            thi = find(fsplist > (levels(thc) / (log(nlist) + 0.5772)));
            if isempty(thi)
                th(thc, 2) = splist(end);
            elseif thi(1) > 1
                th(thc, 2) = splist(thi(1) - 1);
            else
                th(thc, 2) = splist(1) / 2.001;
            end
        end
    end
end
